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Motivated by recent experimental activities on surface critical phe- 
nomena, we present a detailed theoretical study of the near-surface 
behavior of the local order parameter m(z) in Ising-like spin sys- 
tems. Special attention is paid to the crossover regime between 
"ordinary" and "normal" transition in the three-dimensional semi- 
infinite Ising model, where a finite magnetic field H\ is imposed 
on the surface which itself exhibits a reduced tendency to order 
spontaneously. As the theoretical foundation, the spatial behav- 
ior of m(z) is discussed by means of phenomenological scaling ar- 
guments, and a finite-size scaling analysis is performed. Then we 
present Monte Carlo results for m{z) obtained with the Swendsen- 
Wang algorithm. In particular the sharp power-law increase of the 
magnetization, m(z) ~ Hi z 1_,? ° r , predicted for a small Hi by 
previous work of the authors, is corroborated by the numerical re- 
sults. The relevance of these findings for experiments on critical 
adsorption in systems where a small effective surface field occurs is 
pointed out. 

Key words: Surface critical phenomena, critical adsorption, Ising 
model, Monte Carlo simulation 



1 Introduction 



A great deal of current experimental activity concentrates on the investiga- 
tion of critical phenomena near surfaces. After the impressive confirmation 
of the theoretical predictions [1-4] in experiments with binary alloys [5], the 
more recent experimental efforts focus on binary mixtures near their conso- 
lute point [6-8] and near-critical fluids [9]. In most of these experiments the 
order parameter, the concentration difference in fluid mixtures or the density 
difference between liquid and gaseous phase in fluids, plays a central role. 
For instance the reflectivity and ellipticity measured in light-scattering ex- 
periments are directly related to the order-parameter profile [10,11]. Hence, 
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quantitative information about the local order parameter near the boundary- 
is necessary for the interpretation of the experimental data. 

While a very well-developed theory exists for the individual surface univer- 
sality classes, corresponding to fixed point of the renormalization-group flow, 
the picture in the crossover regions between the fixed points is less complete. 
The experiments are generically not carried at the fixed points, however, and 
so a detailed understanding of the crossover region is particularly important. 

Consider for example the semi-infinite three-dimensional (3-d) Ising system 
with spin-spin interaction J. In this model the influence of the surface is 
usually taken into account by means of a modified exchange interaction J\ in 
the surface and a magnetic field H 1 imposed on the surface spins [1,2]. While 
the former models modifications due to the surface within the critical medium, 
the latter represents the influence of the adjacent (noncritical) medium, as the 
container wall for example, on the system. 

For a brief (and necessarily incomplete) summary on surface critical phenom- 
ena, let us first set H\ = 0. Then at the bulk critical point T c the tendency to 
order near the surface can be reduced, increased, or unchanged compared with 
the bulk. Which case is realized depends on the ratio J\/J. At a particular 
value, J{ p ~ 1.5 J [1,12,13] the third case is realized. This corresponds to the 
surface universality class of the "special transition". For J < J{ p the surface 
has a reduced tendency to order but nevertheless becomes (passively) ordered 
at the bulk phase transition. In the opposite case, J > J^ p , the surface orders 
at a temperature above T c , and at T c the bulk undergoes a phase transition in 
the presence of an already ordered surface. From the viewpoint of the renor- 
malization group the special transition is an unstable fixed point [2]. For a 
start value J\ < Jl P the (running) surface coupling is driven to the stable 
fixed point J\ — corresponding to the universality class of the "ordinary 
transition". For J\ > J[ p it is driven to J\ = oo, again a stable fixed point, 
corresponding to the universality class of the "extraordinary transition" [2,14]. 

Next we consider the effects of Hi in a system with J x < J{ p . This is, for 
example, the situation generically met in experiments with binary fluids. In 
particular we are interested in the behavior of the order parameter in this situ- 
ation. The universality classes are determined by the fixed-point values Hi — 
and Hi — oo of the renormalization-group transformations. For H x = 0, at the 
ordinary transition, the order parameter m{z) simply vanishes since, in terms 
of Ising spins, the symmetry under the reversal Sj — > — Sj is not broken, neither 
in the bulk nor in the surface. For Hi = oo the universality class is called the 
"normal transition". The normal transition is known to be equivalent to the 
extraordinary transition [15,16]. In both cases m[z) starts from a large mi 
at the surface and then decays to the bulk equilibrium value (being zero for 
T > T c and nonzero for T < T c ). At T c , i.e. for infinite correlation length £, 
the decay is described by a universal power law m ~ z~Pl v for macroscopic 
distances z. For instance for the 3-d Ising model /3/z/ ~ 0.52 [17]. For T ^ T c 
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a crossover to the exponential decay ~ exp(— z/£) takes place in a distance 
z ~ £ from the surface. 

What happens in the crossover region between ifi = and Hi = oo? Mean- 
field theory predicts a profile that starts from some finite m-y and then monoto- 
nously decays to the equilibrium value. In Ref. [18] the present authors have 
shown that, contrary to the naive (mean- field) expectation, fluctuations may 
cause the order parameter to steeply increase to values m(z) > mi in a 
surface-near regime. This growth is described by a universal power law 

m{z)^H 1 z K with K=l-n°[ d , (1) 

where ?7° rd is the anomalous dimension pertaining to the ordinary transition 
(governing the decay of correlations in the direction perpendicular to the sur- 
face [2]). For instance for the 3-d Ising model k ~ 0.21. 

The scenario for the crossover between ordinary and normal transition devel- 
oped in [18] is the following: At bulk criticality and J\ < J[ p for any finite H\ 
the order-parameter profile increases up to a certain length scale l ord and then 
crosses over to the power law ~ z~Pl v ' . The scale l ord is given by an inverse 
power of Hi and, thus, becomes smaller for increasing H x such that in the 
limit H\ — > oo the maximum has moved to the surface. In this limit only 
the previously mentioned monotonous power-law decay characteristic for the 
normal transition is left over. A qualitative sketch of typical crossover pro- 
files is shown in Fig. 1. In this plot the axes are logarithmic and both m(z) 
and z are measured in arbitrary units. The individual curves have the correct 
asymptotics, m{z) ~ z ' 21 for z — > and m(z) ~ z ~ - 52 for z — > oo. However, 
the (yet unknown) real crossover function is replaced by a simple substitute. 



In Ref. [19] it was demonstrated by means of MC simulations and the compar- 
ison with exact results that also in the 2-d Ising model the crossover between 
ordinary and normal transition is qualitatively of the same form as in d = 3. 
However, the simple power law (1) is modified by a logarithm in d — 2. The 
main purpose of the present work is to verify the results for the order param- 
eter obtained in [18], where scaling and heuristic arguments were used and a 
quantitative calculation in the framework of renormalization-group improved 
perturbation theory was performed, by means of Monte Carlo (MC) simulation 
explicitly for the three-dimensional system. 

The MC studies devoted to (static) critical phenomena near surfaces that are 
contained in the literature concentrated mainly on the dependence of thermo- 
dynamic variables on J\j J and, in particular, on critical adsorption for large 
J\ [20] as well as on the precise location of Jl P [13,21]. H\ was set to zero in 
these studies. Numerical studies of the influence of H\ concentrated on the 
surface layer magnetization [22] and, in the context of wetting, on systems 
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below the critical temperature [23]. To our knowledge, there is no work in the 
literature where order-parameter profiles at or above T c with non-vanishing Hi 
were studied by MC methods and which could have been directly compared 
with the analytic results reported in Ref. [18]. 

The rest of this paper is organized as follows: In Sec. 2 we summarize and 
supplement the main results of [18], especially the phenomenological scal- 
ing analysis which allows to make quite precise predictions for m(z) in the 
crossover regime. In Sec. 3 our MC procedure, essentially the Swendsen-Wang 
algorithm slightly modified to allow for the inclusion of a surface field, is de- 
scribed. The MC results are presented in Sec. 3.3. Eventually, the last section 
contains besides a short summary remarks on the relevance of our results for 
experiments. 



2 Theory 

2.1 Model 



We consider the semi-infinite Ising system with a free boundary on a plane 
square lattice. The exchange coupling between direct neighbors in the bulk is 
J. In the surface the nearest-neighbor coupling is J\. A surface magnetic field 
Hi is imposed on the boundary spins and bulk magnetic fields are set to zero 




Fig. 1. Qualitative sketch of order-parameter profiles in the crossover regime between 
ordinary and normal transtion in double-logarithmic representation. Both z and m 
are given in arbitrary units. 
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such that the Hamiltonian of the model reads 

Rising = - J S i S 3 " J l Y S i S 3 ~ S 'i ' ( 2 ) 

<ij>eV <ij>edV i£dV 

where dV and V stand for the boundary and the rest of the system (without 
the boundary), respectively Below we mainly work with the dimensionless 
variables 

K = J/k B T , K x = J x jk B T , and h x = Hjk B T . (3) 

For the (reduced) critical bulk coupling we took K c = J/k B T c = 0.22165 from 
the literature [17]. The value of K\ that corresponds to the special transition 
is K[ p = 1.5 K c [21,13]. 

2.2 Scaling Analysis 

In the critical regime thermodynamic quantities are described by homogeneous 
functions of the scaling fields. As a consequence, the behavior of the local 
magnetization under rescaling of distances should be described by 

m(z, r, h x ) ~ b~ x * m(zb-\ rb l '\ h x b y ° rd ) , (4) 

where = (3 /is and yl rd = A° rd / v are the scaling dimensions of the equilib- 
rium magnetization m(z — > oo) and the surface field hi, respectively [2]. In 
general the surface exponents have different values for different surface uni- 
versality scales [2], and so these quantities are additionally marked by 1 ord ' 
for belonging to the ordinary transition. The (MC) literature values for the 
3-d Ising model are = 0.518(7) [17] and y° rd = 0.73. The value y° rd was 
obtained by employing the scaling relation yi + x x — d — 1 together with the 
recent Monte Carlo result x°{ d = (3? d /v = 1.27 [24]. 

Removing the arbitrary rescaling parameter b in Eq. (4) by setting it ~ z, one 
obtains the scaling form of the magnetization 

m{z, r, hi) ~ z~ x *M{zH, z/r d ) , (5) 

where 

r d „ h; 1/y ° rd (6) 

is the length scale determined by the surface field. The second length scale 
pertinent to the semi-infinite system and occurring in (5) is the bulk corre- 
lation length £ = t~ v . Regarding the interpretation of MC data, which are 



5 



normally obtained from finite lattices, one has to take into account a third 
length scale, the characteristic dimension L of the system, and a finite-size 
scaling analysis has to be performed. The latter will be described in Sec. 2.6. 

Going back to the semi-infinite case and setting r = 0, the only remaining 
length scale is l ord , and the order-parameter profile can be written in the 
critical-point scaling form 

m(z, h) ~ z~ x * M c {z/l ord ) . (7) 

As said above, for z — > oo the magnetization decays as ~ z~ x <t> and, thus, 
Ai c (() should approach a constant for ( — > oo. In order to work out the short- 
distance behavior of the scaling function Ai c ((), we demand that m(z) ~ m\ 
as z — > 0. This means that in general, in terms of macroscopic quantities, the 
boundary value of m(z) is not mi. If the ^-dependence of m(z) is described 
by a power law, it cannot approach any value different from zero or infin- 
ity as z goes to zero. However, the somewhat weaker relation symbolized by 
"~" should hold, stating that the respective quantity asymptotically (up to 
constants) "behaves as" or "is proportional to" . This is in accord with and ac- 
tually motivated by the field-theoretic short-distance expansion [25,2], where 
operators near a boundary are represented in terms of boundary operators 
multiplied by c-number functions. 

In the case of the 3-d Ising model the foregoing discussion leads to the conclu- 
sion that m(z) ~ hi because the "ordinary" surface with Ki < is param- 
agnetic and responds linearly to a small magnetic field [15]. This is different 
in d — 2 where an additional logarithmic factor occurs [26], mi ~ hi In hi, 
and this logarithm also leaves its fingerprint on the near-surface behavior of 
the magnetization [19]. 

The immediate consequence of the simple linear response in d — 3 for the 
scaling function Ai c {C) occuring in (7) is that it has to behave as ~ ( y ° r in the 
small-C limit. After inserting this in (7), we obtain that the exponent governing 
the short-distance behavior of m(z) is given by the difference between y{ rd and 
#0, such that for z <C l ord the magnetization is described by 

m(z) ~ h z y °i rd - x * . (8) 

Using the scaling relations i]± = (rj + r)\\)/2 and yi = (d — r)\\)/2 [2], we 
eventually obtain k as expressed in (1). 

In the mean-field approximation the result for k is zero. Thus, in this case one 
really has m(z — > 0) = mi and the monotonously decaying order-parameter 
profile mentioned earlier. However, a positive value is obtained when fluctu- 
ations are taken into account below the upper critical dimensionality d* = 4. 
Taking y° rd ~ 0.73 from above, the value for k is 0.21. 
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Phenomena to some extent analogous to the ones discussed above were re- 
ported for the crossover between special and normal transition [27]. Also near 
the special transition the surface field hi gives rise to a length scale. How- 
ever, the respective exponent, the analogy to k, is negative, and, thus, one 
encounters a profile that monotonously decays for all (macroscopic) z, with 
different power laws in the short- distance and the long-distance regime and 
a crossover at distances comparable to the length scale set by hi. However, 
non-monotonous behavior in the crossover region as described above for m(z) 
is a common feature in the case of the energy density in d — 2 [28] as well as 
in higher dimensionality [29]. 

2. 3 Relation to Critical Dynamics 

The spatial variation of the magnetization discussed so far strongly resembles 
the time dependence of the order parameter in relaxational processes at the 
critical point. If a system with nonconserved order parameter (model A) is 
quenched from a high-temperature initial state to the critical point, with a 
small initial magnetization m^, the order parameter behaves asm~ t e 
[30], where the short-time exponent 9 is governed by the difference between 
the scaling dimensions of initial and equilibrium magnetization divided by the 
dynamic (equilibrium) exponent [31]. Like the exponent k in (1), the exponent 
9 vanishes in MF theory, but becomes positive below d*. For example, its value 
in the 3-d Ising model with Glauber dynamics is 9 = 0.108 [32]. 

The high-temperature initial state of the relaxational process is to some ex- 
tent analogous to the surface that strongly disfavors the order and that (for 
hi = 0) belongs to the universality class of the ordinary transition. Further 
expanding this analogy, heating a system from a low-temperature (ordered) 
initial state to the critical point would be similar to the situation at the ex- 
traordinary transition. Eventually, analogous to the special transition would 
be a "relaxation process" that starts from an equilibrium state at T c . 

2.4 Heuristic Argument 

There is also a heuristic argument for the growth of the magnetization in the 
near-surface regime expressed in (1). As said above, a small hi generates a 
surface magnetization mi ~ hi. Regions that are close to the surface will 
respond to this surface magnetization by ordering as well. How strong this 
influence is depends on two factors. 

First, it is proportional to the correlated area in a plane parallel to the surface 
at a distance z. While correlations in the surface are asymptotically (for Ji — > 
0) suppressed, for z > the range of correlations between spins located in a 
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plane parallel to the surface in a distance p from each other can be regarded as 
finite, because for p > z the parallel correlation function is governed by surface 
exponents and decays much faster than in the bulk. Hence the corresponding 
effective correlation length, £\\(z), should behave as z. Referring once more to 
critical dynamics as discussed in the previous section, £\\(z) is analogous to the 
time dependent (growing) correlation length £(£) ~ t 1 ^ (where ( here stands 
for the dynamic equilibrium exponent). 

Second, it depends on the probability that a given spin orientation has "sur- 
vived" in a distance z. For small hi and z < l ord , the latter is governed by the 
perpendicular correlation function C(z) ~ z H . Taking into account 

both factors, we obtain 

m(z) ~ h C(z) ef- 1 = hiz 1 '^, (9) 

the short-distance power law reported in (1). Qualitatively speaking, the sur- 
face when carrying a small mi induces a much larger magnetization in the 
adjacent layers, which are much more susceptible and capable of responding 
with a magnetization m ^> m\. 

This simple picture for the anomalous short-distance behavior holds for dimen- 
sions 2 < d < 4. At and above the the upper critical dimension d* = 4, where 
the mean-field theory starts to provide the correct description, the power-law 
growth of magnetization is not observed, since there the increase of the corre- 
lated surface area is compensated by the decay of the perpendicular correlation 
function. In the case of the two-dimensional Ising model the assumption that 
mi ~ hi is no longer valid, and logarithmic terms occur [19]. 

2. 5 Modifications at T ^ T c 

The phenomenological scaling analysis presented above can be straightfor- 
wardly extended to the case r > 0. In d > 2, we may assume that the behav- 
ior near the surface for z << £ is unchanged compared to (1), and, thus, the 
increasing profiles are also expected slightly above the critical temperature. 
The behavior farther away from the surface depends on the ratio l ord /^. In the 
case of l ord > £ a crossover to an exponential decay will take place for z ~ £ 
and the regime of nonlinear decay does not occur. For l ord < £ a crossover to 
the power-law decay ~ z~^l v takes place and finally at z ~ £ the exponential 
behavior sets in. 

An interesting phenomenon can be observed in the case £ < l ord . As discussed 
above, m(z) then never reaches the regime with power- law decay, but crosses 
over from the near-surface increase directly to the exponential decay. Since 
the region where m(z) grows extends up to the distance £, the magnetization 
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in the maximum has roughly the value m max ~ £ K . Now, the amplitude of the 
exponential decay should behave as ~ m max such that for z ^> £ we have 

m(z)~/iirexp(-z/0. (10) 

In other words, in the case £ < l ord the short-distance exponent k not only 
governs the behavior of m{z) near the surface, but also leaves its fingerprint 
much farther away from the surface in form of an universal dependence of the 
amplitude of the exponential decay on the correlation length ~ £ K . Nothing 
comparable occurs when £ = oo (compare Sec. 2.2 above). When l ord is the only 
scale, all profiles approach the same curve m(z) « Az~^l v for z>l ord , with 
an amplitude A independent of hi. An analogous phenomenon, termed "long- 
time memory" of the initial condition, does also occur in critical dynamics for 
T >T C [33]. 

Below the critical temperature (and near the ordinary transition), the short- 
distance behavior of the order parameter is also described by a power law, 
this time governed by a different exponent, however [34]. The essential point 
is that below T c the surface orders spontaneously even for hi = 0. Hence, 
in the scaling analysis the scaling dimension of h\ has to be replaced by the 
scaling dimension of mi, the conjugate density to hi, given by x°{ d = I3{ rd /v 
[2]. The exponent that describes the increase of the profile is thus x° rd — x^ 
[34], a number that even in mean-field theory is different from zero (= 1) and 
for the 3 — d Ising model its value is 0.75. 

2.6 Finite Size Scaling 

In order to assess the finite size effects to be expected in the MC simula- 
tions, we have to take into account the finite-size length scale L. The latter 
is proportional to the linear extension of the lattice (called N below). The 
generalization of (4) reads [35] 

m(z, r, hi,L) ~ b~ x * m{zb~\ rb 1/u , hi b v ° r \ Lb' 1 ) , (11) 

and proceeding as before, we obtain as the generalization of (5) to a system 
of finite size: 

m(z, r, hi, L) ~ z- x +M(z/£, z/l ord , z/L) . (12) 

Thus even at T c there are two macroscopic length scales, on the one hand L 
(imposed by the geometry) and on the other hand l ord (the scale set by hi). 

It is well known that for large z >L we have to expect an exponential decay of 
m(z) on the scale L. In the opposite limit, when z is smaller than both L and 
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l ord , we expect the short- distance behavior (1) to occur. That this expectation 
turns out to be correct is the necessary condition for observing (1) in MC 
simulations. Stated in terms of correlation lengths it means that as long as £y 
(cf. discussion in Sec. 2.4) is smaller than L (and the bulk correlation length 
£), the form of the profile is unchanged compared to the one of the semi- 
infinite system (at bulk criticality). In particular it implies that the surface 
magnetization mi, whose linear response to hi was an important ingredient 
to our scaling analysis of Sec. 2.2, should not depend on L, as long as L can 
be regarded as macroscopic. 

Farther away from the surface, the form of the profile depends on the ratio 
between l ord and L. In the case of l ord > L a crossover to an exponential decay 
will take place for z ~ L, and, analogous to the situation with a finite corre- 
lation length, also in the finite-size system the amplitude of this exponential 
decay is governed by the exponent k (compare to (10) and the discussion in 
Sec. 2.3), such that we have for z 3> L 

m(z) ~ h L K exp(-z/L) (13) 

Again, analogous finite-size effects were reported also in relaxation processes 
near criticality [33]. In the opposite case, l ord < L, a crossover to the power- 
law decay ~ z~ x * takes place, followed by the crossover to the exponential 
behavior at z ~ L. Thus, qualitatively, the discussion for systems of finite size 
is largely analogous to the one in Sec. 2.3 for finite £. 



3 Monte Carlo Simulation 

3. 1 System 

The results of the scaling analysis, especially the short-distance law (1), were 
checked by MC simulations. To this end we calculated order-parameter profiles 
for the 3-d Ising model with uniform bulk exchange coupling K and set K\ = 0, 
corresponding to the fixed-point value of the ordinary transition. 

The geometry of the systems studied was that of a rectangular (cuboidal) 
lattice with two free surfaces opposite to each other and the other boundaries 
periodically coupled. The surface field hi was imposed on both free surfaces. 
The linear dimension perpendicular to the surfaces was taken to be two times 
larger than the lateral extension in order to keep corrections due to the second 
surface, the so-called Fisher-de Gennes effect [36], small [36]. Hence, when we 
talk about a lattice of size N in this section, we refer to a rectangular system 
with N 2 x 2N spins. 
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The distance from the surface is still called z in the following, although it 
is clearly an integer quantity, with z = corresponding to the location of 
one of the surfaces. Order parameter profiles were calculated by averaging in 
individual configurations over planes parallel to the surface and, in turn, we 
averaged over a large number of configurations generated by the algorithm 
described in Sec. 3.2. Eventually the symmetry of the system was used and 
also the average between the left and right half of the lattice was taken. 

3. 2 Procedure 

We consider the Ising model, defined by (2). The aim of the (equilibrium) 
Monte Carlo procedure is to generate a representative sample of configurations 
s distributed according to the Boltzmann factor P(s) ~ exp [H(s)/ksT] [37]. 
Further, it must be guaranteed that, starting from any initial configurations, 
after a reasonable amount of time such a sample can be extracted. The latter 
is in principle provided if the algorithm that generates a new configuration s' 
from the old one satisfies detailed balance. In terms of transition probabilities 
W(s — >s') this condition can be expressed as 

W(s->s') exp [-H(s)/k B T] = W(s'-> s) exp [-H(s')/k B T] . (14) 

For practical purposes, however, not any algorithm satisfying (14) is suitable 
for MC simulations of critical or near-critical systems. The reason is that phys- 
ically meaningful algorithms, like Glauber and Kawasaki dynamics [37], are 
greatly hampered by the critical slowing down upon approaching the equilib- 
rium. One way out of this dilemma is the Swendsen-Wang (SW) algorithm 
[38], which satisfies (14) but does (probably) not correspond to a physically 
meaningful dynamics. 

The SW algorithm generates a transition (or update) s — > s' between spin 
configurations via connected bond clusters. A cluster configuration n is con- 
structed from s by creating bonds between neighboring spins of equal sign. 
Then these bonds are "activated" [39] with probability 

p=l-e- 2K . (15) 

No bonds are generated between spins of opposite sign. As the next step, bond 
clusters are defined as connected sets of active bonds. Also isolated spins are 
identified as a cluster, such that eventually each spin belongs to a cluster. 

In order to obtain a new spin configuration s' from n, one assigns to all sites of 
a given cluster a new spin value with equal probability for each spin direction 
(independent of the old spin value). The probability for the transition s — > s' 

W(s|n|s') =p\l -p) m q- Nc (16) 
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where n is an intermediate cluster configuration with N c clusters, and b and m 
are the numbers of "active" and "inactive" bonds, respectively. This transition 
corresponds to one Monte Carlo sweep. 

In order to verify that the algorithm satisfies detailed balance, we have to 
consider a transition in the opposite direction. For the transition s' from s via 
the same cluster configuration n, there is a probability 

W(s|n|s') = p\l - p) m 'q- Nc , (17) 



with the same b and N c as before. However, the number of non-active bonds 
m' can in general be different, because neighboring clusters can originate from 
domains with spins of the same or of different sign, in both cases leading to 
the same cluster configuration. 

The total transition probability from s to s' is given by 

W( s ,s') =^W(s|n|s'), (18) 



where the sum runs over all possible intermediate cluster configurations n. 
Since the sum b + m is constant for a given spin configuration, it is straight- 
forward to show that 

Eventually, taking into account that the energy difference between s and s' is 
given by 

- 2 J(m - m!) = -AH , (20) 



with (15) one obtains the detailed balance relation (14). 

The algorithm presented so far works as long as no magnetic fields are imposed 
on the spins. To take into account the third term in (2) that describes the 
influence of the surface magnetic field Hi, we follow Wang [40] and introduce 
a layer of "ghost" spins next to the surface that couple to the surface spins 
only. The "ghost" spins all point in the direction of Hi and couple to the 
"real" spins with coupling strength equal to Hi. If one or more "active" bond 
between a surface and a ghost spin exist, the cluster has to keep its old spin 
when the system is updated. This prescription preserves detailed balance. In 
the practical calculation this rule was realized by a modified (reduced) spin-flip 
probability 

p(n s ) = 1 - i exp(-2 hi n s ) (21) 
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for clusters pointing in the direction of hi, where n s is the number of surface 
spins contained in the cluster. For clusters pointing in opposite direction the 
probability has to remain unchanged (equal to 1/2). 

In order to obtain an equilibrium sample, we discarded several hundred - 
the precise number depended on the size of the system - configurations after 
the start of the run. To keep memory consumption low, we used multispin- 
coding techniques, i.e. groups of 64 spins were coded in one long integer. All 
calculations were run on a Silicon Graphics computer (Power Challenge) with 
four Risk 8000 processors. To obtain a profile with reasonable statistics for 
our largest system (N=256) took about one week of (single-processor) CPU 
time. 



3.3 Results 



From the magnetization profiles especially the surface magnetization mi as a 
function of hi can be extracted. It is instructive to compare the results for the 
3-d Ising model with those for the two-dimensional case obtained in Ref. [19]. 
This is done in Fig. 2. The crosses represent the data obtained from a three- 
dimensional system with N = 256, the circles stem from the two-dimensional 
Ising model with lattice size 512x2048. In both cases the statistical errors are 
smaller than the symbol size. 




Fig. 2. Surface magnetization mi as a function of hi for the Ising model in d = 2 
(full circles) and d = 3 (crosses) in double-logarithmic representation. The data for 
d = 2 stem from a 512x2048 lattice, those for d = 3 from a 256 2 x 512 system. The 
solid line shows the exact result for the d = 2 semi-infinite system. The dashed line 
is a fit to the our MC data. 
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The situation in three dimensions is obviously quite simple. Up to hi ~ 0.1 
the response of m\ on hi is just linear. This is the regime where the scaling 
analysis of Sec. 2.2 applies, in particular the basic assumption that m(z) ~ hi 
as z goes to zero. The dashed line is a linear fit to the data for hi < 0.1. For 
larger values of hi the surface magnetization saturates, such that for hi — > oo 
the curve has to approach unity 

The dependence of mi on the surface field in two dimensions is more com- 
plicated. As discussed in detail in Ref. [19] and investigated in many exact 
calculations [26], there occurs a logarithmic factor in the functional depen- 
dence of mi on hi. The solid line shows the exact result of the semi-infinite 
system, which, due to the logarithm, never goes through a regime of linear be- 
havior. The MC data (see Ref. [19]) deviate from the exact curve for hi <0.02 
and for small hi indeed show a linear dependence. This observation, linear 
response for small hi and an approach to the true semi-infinite behavior for 
larger hi, is a finite-size effect consistent with exact results [26]. 

The next point are finite-size effects in the vicinity of the surface, especially 
concerning the dependence of mi on N. As discussed in Sec. 2.6, we expect 
that mi and the profile up to a certain distance ~ iV should not depend on 
N. In Fig. 3 we plotted the data for the local magnetization for four different 
system sizes ranging from N = 64 to 256 up to z = 20. In all cases the surface 
field was hi = 0.01. Quite obviously, mi itself does not vary with N in the 
given range of sizes, confirming the finite-size scaling analysis in Sec. 2.6 as 
well as the assumptions underlying the scaling analysis in Sec. 2.2. From the 
given value of mi all profiles increase for z increasing away from the surface. 
For the first few layers the curves lie on top of each other, but in the smaller 
systems the slopes become smaller compared to larger N at relatively small 
distances already. For the system with 64 2 x 128 spins the regime with growing 
magnetization extends to z ~ 7 only. 

0.028 
0.026 

m (z) 




10 12 14 16 18 20 

Z 



Fig. 3. The local magnetization at and in the immediate vicinity of the surface for 
hi = 0.005 and N = 64 (dashed-dotted), 128 (dotted), 196 (dashed), and 256 (solid 
line) . 
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Fig. 4. Order-parameter profiles for the same parameters as in Fig. 3 in dou- 
ble-logarithmic representation. 

The variation of the same curves on larger scales is displayed in Fig. 4 in 
double-logarithmic form. Going to larger system size the distance of the max- 
imum z m ax grows roughly as ~ N. For N = 256 we have z max ~ 30. Recalling 
the results of the finite-size scaling analysis of Sec. 2.6, we conclude that with 
these parameters the model is in the regime where L < l ord , and going to 
a smaller hi that would increase l ord would not help to extend the region of 
growing magnetization. For z — > the form of the profiles is consistent with 
a power law. However, in the small systems the finite size effects cause the 
profiles to crossover to the exponential decay (compare Sec. 2.6) at a rather 
small distance. Even in the largest system N = 256 that could be studied 
by our present means with reasonable effort, the near-surface power law does 
not extend beyond z ~ 20. The problem is that we indeed have L ~ N, but 
apparently with a rather small constant of proportionality. 

Nevertheless a rough value for the short- distance exponent can be extracted 
from the profile for N = 256. The result is k = 0.16(2). This is somewhat lower 
than our expectation. The deviation from the expected value 0.21 is very likely 
due to the finite-size effects. We can not claim to see the power law (1) over a 
really macroscopic range before the crossover to the finite-size (exponential) 
behavior sets in. So the determination of a more reliable value of k from the 
short-distance behavior remains as a task for larger-scale simulations. 

Magnetization profiles for iV = 256 and hi varying in a wide range between 
0.005 and 5 are plotted in Fig. 5. The dashed line represents the pure power law 
~ £-0.52 cnarac teristic for the extraordinary or normal transition. For small 
hi, up to hi ~ 0.1, the curves show the near-surface growth consistent with 
(1). For hi up to about 0.02, the location of the maximum (here z max ~ 30) is 
determined by the finite-size scale L. Setting hi to larger values, the maximum 
moves closer to the surface. This is the regime where l ord is smaller than L 
and the location of the maximum is governed by l ord . In the case of hi >0.2, 
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the profiles decay monotonously. Setting hi = 5.0, the magnetization at the 
surface is very close to one, and the decay for 10<2;<100 is consistent with 
the power law ~ z~ x< t> The value of the exponent = (3/v obtained from this 
curve is 0.51(1), which is in excellent agreement with the literature value 0.517 
[17]. The up-bending of the profiles for z >100 is due to the second surface. 

Eventually, Fig. 6 shows the values of z max as a function of hi determined 
from the profiles of Fig. 5. From these data, the qualitative picture from above 
can be made more quantitative. Especially with the result Eq. (6) we can in 
principle determine directly the scaling dimension y° rd . A power law fitted 
to the data for 0.02 < hi < 0.1 yields l/yl rd = 1-4(1) (where the error was 



m (z) 



0.1 



0.01 




Fig. 5. Order-parameter profiles for N = 256 and hi = 0.005, 0.01, 0.02, 0.04, 0.08, 
0.1, 0.2, and 5.0 in double-logarithmic representation. The pure power law ~ z -o.52 
(dashed line) is shown for comparison. 




Fig. 6. z max in dependence of h\ in double-logarithmic representation as obtained 
from the profiles of Fig. 5. The error bars were estimated. The dashed line depicts 
the derived power-law dependence of the length scale l ord on hi. 
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estimated), which, in turn, yields 0.71(4) for the scaling dimension y° rd . As 
mentioned above the literature value is 0.73. However, as in the case of the 
short- distance exponent, we also here have to admit that we are not really in a 
regime where we can call l ord large compared with the lattice spacing, and the 
good agreement with what we expected from the scaling analysis is actually 
surprising. 



4 Summary and Concluding Remarks 

We studied the near-surface behavior of the order parameter in the three- 
dimensional Ising system under the influence of a surface magnetic field H ± . 
The anomalous behavior found in [18] by employing scaling arguments and 
perturbative methods was confirmed in the present work by Monte Carlo sim- 
ulations. Especially, the short-distance power law (1) was corroborated. 

However, in our Monte Carlo study especially the region with small Hi is 
severely affected by finite size effects. Even in our largest system (256 2 x 512) 
the increase of m(z) does not extend beyond z ~ 30. In order to obtain reliable 
results for the exponent k determined with the help of the short- distance 
power law (1) and profiles that can be used for the quantitative comparison 
with experimental data, one has to go to systems beyond the size that we are 
able to treat by our present means. 

Concerning experiments on surface critical phenomena, our results should be 
of interest especially in those cases where a small Hi occurs, at a surface that 
disfavors the order. An example where this was obviously realized is the sys- 
tem studied by Desai et al. [8]. In their experiment a binary fluid was studied 
in a container whose walls as a function of time change their preference from 
one component to the other, the time scale of this change being of the order 
of days. In other words, the surface field Hi changes sign as time goes by, and 
during a certain period Hi is small. First substantial steps towards a theo- 
retical explanation of this and other similar experiments on binary mixtures 
were made already by Ciach et al. [41]. A complete theoretical analysis would 
require a careful derivation of experimentally observably quantities like reflec- 
tivity and ellipticity for light scattering experiments on the basis of our results 
for the order-parameter profile. 

Another experiment, discussed already in some detail in [18], is the one by 
Mailander et al. [5] on Fe 3 Al. This system undergoes (among other transi- 
tions) a second-order phase transition between a phase with D0 3 structure 
and one with B2 structure. The near-surface regime was studied by scattering 
of evanescent x-rays. The exponents observed were consistent with the ex- 
pectation for the ordinary transition, but Bragg peaks revealed the existence 
of long-range order near the surface reminiscent to the normal transition. In 
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order to explain the experimental results of [5] on the basis of our findings we 
have to assume that there exits an effective Hi in this system. Then, if the 
associated length scale l ord is larger than both the scattering depth and the 
bulk correlation length, the structure function is governed by the anomalous 
dimension of the ordinary transition. On the other hand, the steep increase 
of the order parameter should provide the explanation for the observed long- 
range order near the surface. 
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